#include "qgaus.h"
#include <math.h>

//--------------------------------------------------------------------
// After debugging/testing, you can put the template stuff in
// its own header

//--------------------------------------------------------------------
// Start of template stuff
// Derived from the quad3d_t in post number 3 of this thread:
//http://www.nr.com/forum/showthread.php?t=1663
//
template < class T >
struct NRf2
{
    Double xsav;
    T *func2d;
    Double operator() (const Double y) {
        return (*func2d) (xsav, y);
    }
};

template < class T, class Y1, class Y2 >
struct NRf1 {
    Y1 & y1;
    Y2 & y2;
    NRf2 < T > f2;
    NRf1(Y1 & yy1, Y2 & yy2):y1(yy1), y2(yy2) { }
    Double operator() (const Double x) {
        f2.xsav = x;
        return qgaus(f2, y1(x), y2(x));
    }
};

template < class T, class Y1, class Y2 >
Double quad2d_t(T & func, const Double x1, const Double x2, Y1 & y1, Y2 & y2)
{
    NRf1 < T, Y1, Y2 > f1(y1, y2);
    f1.f2.func2d = &func;
    return qgaus(f1, x1, x2);
}
// End of template stuff
//--------------------------------------------------------------------

// Test program for 2D integration



struct Ftor
{   double row;
	Ftor(const Double & sphr):row(sphr) {}
    Double operator() (const Double & x, const Double & y) const
	{
	double b;
		if ((row*row - x*x - y*y)<0)
		{
		// ShowMessage(row*row - x*x - y*y);
		 b= sqrt(fabs(row*row - x*x - y*y));
		}
		else
		b= sqrt(row*row - x*x - y*y);


		return b;
    }
};

//
// Use functors to pass functions to quad2d.
//
// This way the radius can be incorporated in the
// calculations without making it a global variable
//
//
// For a given value of x, this calculates
// the lower limit for integration over y.
//
struct Fty1
{
//	double row,r,d;

//	Fty1(const Double & r):radius(r) {}

	Double operator() (const Double & x) const
	{
        return 0;
    }
};


//
//For a given value of x, this calculates
// the upper limit for integration over y.
//

struct Fty2
{
	double row,r,d;

	Fty2(const Double & cylr,const Double & sphr,const Double & dcsphacyl):row(sphr),r(cylr),d(dcsphacyl) {}
	void initialise(const Double & cylr,const Double & sphr,const Double & dcsphacyl)
	 {	row=sphr;
		r=cylr;
		d=dcsphacyl;
	 };
	Double operator() (const Double &x) const
	{
		return GetMin(sqrt(fabs(r*r - (x-d)*(x-d))),sqrt(fabs(row*row - x*x)));
    }
};
/*
//
// Some compiler library math.h headers define M_PI, some do not
//
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif


//
// The integral of r**2 over a circle with a given radius
// is obtained easily with polar coordinates.  Here's
// the analytical result.
//
Double actual(Double r)
{
    return 2.0*M_PI * pow(r,4)/4.0;
}

int main()
{
    Double num_values = 10;
    cout << "Integration of r^2 over a circle using quad2d_t"
         << endl << endl;

    cout << setw(11) << "Radius" << setw(15) << "Approx";
    cout << setw(15) << "Actual" << endl;
    cout << scientific << setprecision(6);

    for (int i = 1; i <= num_values; i++) {

        Double radius = Double(i)/num_values;
        Double xmax   = radius;
        Double xmin   = -xmax;

        // Instantiate functors
        Ftor func;
        Fty1 y1(radius);
        Fty2 y2(radius);

        // Calculate the integral for this radius
        Double result = quad2d_t(func, xmin, xmax, y1, y2);
        cout << setw(15) << xmax
             << setw(15) << result
             << setw(15) << actual(radius) << endl;

    }

    return 0;
}
*/
